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Magnetic monopoles have eluded experimental detection since their 
prediction nearly a century ago by Dirac [1]. Recently it has been shown 
that classical analogues of these enigmatic particles occur as excitations 
out of the topological ground state of a model magnetic system, dipolar 
spin ice ^2 1 _ These quasi-particle excitations do not require a modification 
of Maxwell's equations, but they do interact via Coulombs law and are of 
magnetic origin. In this paper we present an experimentally measurable 
signature of monopole dynamics and show that magnetic relaxation measu- 
rements in the spin ice material Dy^Ti^O-j [3] can be interpreted entirely in 
terms of the diffusive motion of monopoles in the grand canonical ensemble, 
constrained by a network of "Dirac strings" filling the quasi-particle vacuum. 
In a magnetic field the topology of the network prevents charge flow in the 
steady state, but there is a monopole density gradient near the surface of 
an open system. 



Spin ice systems (4 (| such as Dy 2 Ti 2 0j and Ho2Ti 2 Of, can be described by a corner 
sharing network of tetrahedra forming a pyrochlore lattice of localized magnetic mo- 
ments, as shown in figure The pairwise interaction is made up of both exchange and 
dipolar terms 
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where the rare earth ions carry a moment of ten Bohr magnetons, m ~ and 
where Sj is a spin of unit length. The coupling constants are on the 1 K energy scale ; 
for example for Dy^Ti^^i |J|m 2 pa 3.72 K and Dm 2 pa 1.41 if [7J. These energy scales 
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are 100 times smaller than the crystal field terms |8| that confine the spins along the 
axis joining the centres of two adjoining tetrahedra. As a result, on the 1 K energy 
scale the moments behave as Ising spins along this axis. Remarkably, within this Ising 
description the long ranged dipolar interactions are almost perfectly screened at 
low temperature, with the result that the low energy properties are almost identical 
to those of an effective frustrated nearest neighbour model with antiferromagnetic 



interactions 



lOj of strength J e ff = (5D — J)m 2 /3 . This is equivalent to Pauling's model 



for proton disorder in the cubic phase of ice which has extensive ground state 
entropy and violates the third law of thermodynamics |6j. It successfully reproduces 
the thermodynamic behaviour of both ice jl2j and spin ice [13] and describes the 
microscopic properties of the latter to a good approximation. The extensive set of 
spin ice states satisfy the Bernal - Fowler ice rules 14j ; a 3d analogue of the 6 vertex 
model with topological constraint consisting of two spins pointing into and two out of 
each tetrahedron (2 in - 2 out), as shown in figure [1^1. Flipping one spin breaks the 
constraint leaving neighbouring tetrahedra with 3 in - 1 out and with 3 out - 1 in, which 
constitute a pair of topological defects. Within the nearest neighbour model, creation of 
the defect pair costs energy 4 J e ff, while further spin flips can move the defects at zero 
energy cost. It has recently been shown {2} that including the full dipolar Hamiltonian 
of equation (1) leads to an effective Coulombic interaction between the topological 
defects separated by distance r , /io'Mj / 'Airr , where /io is the permeability of free space, 
qi = ±q = ±2m/a, and a is the distance between two vertices of the diamond lattice (see 
figure [I]) ; that is, to a Coulomb gas of magnetic monopoles. Standard electromagnetic 
theory does allow for such excitations |l5|, which correspond to divergences in the 
magnetic intensity H, or magnetic moment M, rather than in the magnetic induction : 
V • B = V ■ (H + M) = 0. On all length scales above the atomic scale, a 3 in - 1 out 
defect appears to be a local sink in the magnetic moment and therefore as a source of 
field lines in H . It can lower its energy by mov ing in the direction of an external field 



16] . What is remarkable about spin ice 



and therefore carries a positive magnetic charge 
is that it allows for the deconfinement of these effective magnetic charges so that they 
occur in the bulk of the material on all scales, rather than just at the surfaces within 



a coarse grained description 



15|. A two dimensional equivalent may exist in artificial 



spin ice, constituting arrays of nanoscale magnets 
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Given the accessibility of these magnetic quasi-particles, the development of an 
experimental signature is of vital importance and interest. The "Stanford" superconduc- 
ting coil experiment could in principle detect the passage of a single magnetic 
quasi-particle, but this seems highly unlikely given that the charges have no mass and 
therefore have diffusive, rather than Newtonian dynamics. A more promising starting 
point is therefore to look for a monopole signal from magnetic relaxation of a macrosco- 



pic sample 



19] . The general dynamic behaviour of spin ice is illustrated in figure [2] 
by the magnetic relaxation time, as a function of temperature for Dy2Tt20j {3}, taken 
from bulk susceptibility measurements. The energy scales discussed above give rise to 
different regimes : the time scale increases in the thermally activated high temperature 
regime, entering a quasi-plateau region below 12 K associated with quantum tunnelling 
processes 8|, before experiencing a sharp upturn below 2K. The spins are Ising like 
below 12 K and the configuration evolves by quantum tunnelling through the crystal 
field barrier, while above this temperature higher crystal field levels are populated and 
the time scale drops dramatically. The quantum tunnelling plateau regime can therefore 
be well represented by an Ising system with stochastic single spin dynamics and hence 
should be dominated by the creation and propagation of monopole objects. This is 
illustrated, in a first approximation, by comparing the data with an Arrhenius law 
t = tq exp(2 Jeff/ksT), as shown by the red curve in figure [2j The time scale To is fixed 
by fitting to the experimental time at AK with J e ff = 1.11 K, the value estimated for 
Dy2Ti 2 0j 3|. 2 J e ff is the energy cost of a single, free topological defect in the nearest 
neighbour approximation and is half that for a single spin flip. The calculation fits the 
data over the low temperature part of the quasi-plateau region, where one expects a 
significant defect concentration without any double defects (4-in or 4-out), and gives 
surprisingly good qualitative agreement at lower temperature, as the concentration 
decreases. Although still in the tunnelling regime, the plateau region corresponds 
to high temperature for the effective Ising system. Good agreement here provides a 
stringent test and any theory not fitting must be discarded. The above expression 
clearly does a good job, allowing us to equate To with the microscopic tunnelling time. 
This test therefore already provides very strong evidence for the fractionalization of 
magnetic charge Q and the diffusion of unconfined particles. However, this (or any 
other) Arrhenius function ultimately fails, underestimating the time scale at very low 
temperature : while it is possible to fit the data reasonably below 2 K by a single 
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exponential function by varying the barrier height, simultaneous agreement along the 
plateau and at lower temperature is impossible. The role of the missing Coulomb 
interaction is therefore clear : although non-confining it must considerably increase the 
relaxation time scale by modifying the defect concentration and slowing down diffusion 
through the creation of locally bound pairs. 



We have tested this idea by directly simulating a Coulomb gas of magnetically charged 
particles (monopoles), in the grand canonical ensemble, occupying the sites of the dia- 
mond lattice. The magnetic charge is taken as q$ = ±q. In the grand canonical ensemble 
the chemical potential is an independent variable, whose value in the corresponding 
magnetic experiment is unknown. In a first series of simulations we have estimated it 
numerically by calculating the difference between the Coulomb energy gained by crea- 
ting a pair of neighbouring magnetic monopoles and that required to produce a pair 
of topological defects in the dipolar spin ice model, with parameters taken from refe- 
rence J2I, giving a configurationally averaged estimate [i/kg = 8.92 K. In a second series 
of simulations [i was taken as the value required to reproduce the same defect concentra- 
tion as in a simulation of dipolar spin ice at temperature T. Here fi varied by 3% only, 
with the same mean value as in the first series, showing that our procedure is consistent. 
The chemical potential used is thus not a free parameter. As the Coulomb interaction 



is long ranged, we treat a finite system using the Ewald summation method [20l . |21| . 
The monopoles hop between nearest neighbour sites via the Metropolis Monte Carlo 
algorithm, giving diffusive dynamics, but with a further local constraint : in the spin 
model a 3 in - 1 out topological defect can move at low energy cost by flipping one of the 
three in spins, the direction of the out spin being barred by an energy barrier of 8 J e ff- 
An isolated monopole can therefore hop to 3 out of 4 of its nearest neighbour sites only, 
dictated by an oriented network of constrained trajectories similar to the ensemble of 
classical "Dirac string" |2| of overturned dipoles 15|. The positively charged monopoles 
move in one sense along the network while the negative charges move in the opposite di- 
rection (see figure db). The network is dynamically re-arranged through the evolution of 
the monopole configuration. The vacuum for monopoles in spin ice thus has an internal 
structure ; the Dirac strings which, in the absence of monopoles, satisfy the ice rules at 
each vertex. This structure is manifest in the dynamics and influences the resulting time 
scales. In fact the characteristic time scale that we compare with experiment comes from 
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the evolution of the network of Dirac strings rather than from the monopoles themselves. 
Indeed, the monopole autocorrelation time, as extracted from the monopole density - 
density correlation function 22J turns out to be small for this range of temperature. 



We locally define the string network by an integer a = ±1, giving the orientation of 
the Dirac string along each bond of the diamond lattice and define the autocorrelation 
function 

C ® = ^X>(*M°), ( 2 ) 

i 

where t is the Metropolis time and N is the number of bonds (up to iV w 25000). 
For the initial conditions we take an ordered network with no monopoles, which we 
let evolve at temperature T until an equilibrium configuration is attained. This defines 
t = 0. C(t) decays almost exponentially, with characteristic time, r, that varies with 
temperature. To avoid initial transient effects we define r such that C(r) = 0.8. The 
time is re-set to zero when C(t) decays beyond 0.01 and the process is repeated many 
times to give the configurationally averaged decay time. In figure [2] we compare our 
simulations with the experimental data of reference [3|. The Metropolis time is again 
scaled to the experimental time at 4 if and there is again no scale factor on the 
temperature axis. Data for fixed chemical potential are shown by the pink triangles, 
while data with /x varying are shown by the blue circles. There is a quantitative 
evolution of the simulation data compared to the nearest neighbour spin ice model. 
Agreement between the experimental and numerical data now looks excellent, showing 
clearly that the experimental relaxation is due to the creation and proliferation of 
quasi-particle excitations that resemble classical monopoles in the magnetic intensity H. 
As the temperature increases, towards the end of the plateau region, a small systematic 
difference occurs. This is because the spin system can access double defects at finite 
energy cost, while this state corresponds to two like charges superimposed on the same 
site, which is excluded by the Coulomb interaction. The inset of figure [2] shows results 
at low temperature illustrating in detail the extent of the improvement in comparison 
with experiment. Allowing the variation of \i provides a further evolution towards the 
experimental data, compared with that for fixed /i and the blue circles represent our best 
numerical results. We now have quantitative agreement between experiment and theory 
down to low temperature, showing that the Coulomb interactions are responsible for the 
non-Arrhenius temperature dependence of the relaxation time scales. Differences remain 
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at this level of comparison below 1 K but to go further would require an even more 



detailed modelling of spin ice 



23| as well as complementary experimental measurements. 



Finally we consider the response of monopoles to an external magnetic field, h, pla- 
ced along one of the [100] directions. Applying such a field to a system for closed circuit 
geometry (periodic boundaries), one might expect the development of a monopole cur- 
rent in the steady state [if]]. This is not the case, at least for the nearest neighbour 
model, where we find that a transient current decays rapidly to zero (see figure [3^). 
The passage of a positive charge in the direction of the field re-organises the network 
of strings, leaving a wake behind it that can be followed either by a negative charge, or 
by a positive one moving against the field, with the result that the current stops. This 
is a dynamic rather than static effect and is not related to confinement of monopole 
pairs by the background magnetization [2|. Reducing the temperature at finite field, 
the magnetization saturates around a critical temperature ; a vestige of the Kasteleyn 



transition 24] , which is unique to topologically constrained systems. Confinement occurs 
here, as the Zeeman energy outweighs the entropy gain of free monopoles. The transient 
currents suggest the development of charge separation in an open system. This is indeed 
the case despite the fact that monopole numbers are not conserved at open boundaries. 
In figure [3b we show the profile of positive charge density across a sample of size L, with 
open boundaries, for varying fields. There is a clear build up of charge over a band of 4-5 
lattice spacings, although including long range interactions may lead to a quantitative 
change in this value. As the ratio T/h and the monopole density go to zero the band 
narrows and the system forms a conventional layer of magnetic surface charge as one 
expects for any magnetically ordered system [15]. In the absence of topological defects 
the magnetization is conserved from one layer to another, so that a charge density profile 
manifests itself as a magnetization profile. The data here suggest charge build up in a 
layer several nanometres thick, making it in principle a measurable effect. 
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Figure 1: Spin Ice structure and emergence of monopoles. a) The magnetic ions (Ho 3+ 
or Dy 3+ ) lie on the sites of the pyrochlore lattice and are constrained to the bonds of the 
duel diamond lattice (dashed lines). Local topological excitations 3 in - 1 out or 3 out - 1 in 
correspond to magnetic monopoles with positive (blue sphere) or negative (red sphere) charges 
respectively, b) The diamond lattice provides the skeleton for the network of Dirac strings with 
the position of the monopole restricted to the vertices. The orientation of the Dirac strings 
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FIGURE 2: Relaxation time scales r in Dy2Ti207 : experiment and simulation. The 

experimental data (x) are from Snyder & al. [3]. The Arrhenius law (red line) represents the 
free diffusion of topological defects in the nearest neighbour model. The relaxation time scale 
of the Dirac string network driven by Metropolis dynamics of magnetic monopoles has been 
obtained for fixed chemical potential (pink A) and with /i varying slowly to match the defect 
concentration in dipolar spin ice (blue ©). The temperature scale is fixed without any free 
parameters. Inset : Same data shown in the low temperature region. 
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FIGURE 3: Monopole density profile, a) A magnetic field h is applied at t = along the 
z-axis ([100] direction). We display the transient flux of positive charges, <3?, passing through 
a plane perpendicular to the field, as a function of Metropolis time, t. The simulations are 
obtained using the nearest neighbour spin ice model with periodic boundary conditions (H) 
and open boundaries, with current measured either at the surface (blue A) or in the bulk (red 
9). b) Density of positive defects in the horizontal planes for T = IK, as a function of z and 
h (in units of k^T '/mfMo)- 



